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AUTOMATIC  PIECEWISE  LINEARIZATION  IN  IDEAL  PLASTICITY^ 

by 

Philip  G.  Hodge, 


ABSTRACT 


A method  is  proposed  for  constructing  a piecewise  linear  approximation 
to  an  arbitrary  yield  function.  The  approximation  is  constructed  in  the 
course  of  solving  a given  boundary-value  problem.  Its  use  is  illustrated 
with  some  simple  examples. 
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AUTOMATIC  PIECEWISE  LINEARIZATION  IN  IDEAL  PLASTICITY 


1.  Introduction.  The  defining  equations  of  a typical  boundary  value 
problem  for  a linear-elastic/perfectly-plastic  material  consist  of  the 


equilibrium  equations 


a . . . + F . = 0 
,J  1 


the  strain-displacement  relations 


"ij  " 


and  the  constitutive  equations 


S'j  " 


f < f. 


X > 0 


(3b, c) 


IF  (f  < f(3  OR  f < 0)  THEN  A = 0 (3d) 

together  with  sufficient  prescribed  boundary  conditions  on  stress  or 
displacement.  In  Eqs.  (3),  the  yield  function  f is  a given,  closed,  convex 
function  of  the  stresses,  f^  is  a given  constant,  are  the  elastic 

constants,  X is  an  unknown  scalar  variable,  and  a superposed  dot  means 
differentiation  with  respect  to  any  time-like  parameter. 

The  usual  method  of  solving  this  problem  is  to  assume  that  a complete 
solution  is  known  at  a gene>"ic  time  t^,  to  solve  the  rate  problem  at  t^, 
and  then  to  integrate  ahead  to  time  t^  + At.  The  rate  problem  is  defined 
by  the  differentiated  forms  of  (1)  and  (2),  together  with  (3)  which  may  be 
rewritten  as  (3a)  together  with 

IF  f < f^  THEN  A = 0 

• • 

ELSE  IF  f < 0 THEN  A = 0 (4) 

• • 

ELSE  (f  = 0 AND  A > 0) 


2. 


We  observe  that  the  first  branch  in  (4)  is  defined  by  the  known  solution 

t^,  but  that  the  correct  second  branch  is  defined  only  in  terms  of  the 

rate  solution.  Difficulties  associated  with  choosing  the  correct  second 

branch  have  been  discussed  elsewhere  [1,2]  ; here  we 

shall  assume  that  we  have  already  identified  the  correct  branch.  Once  the 

rate  problem  is  solved,  the  stress  at  time  t^  + At  may  be  written  approximately 

as 

‘ “ “ij  (5> 

We  are  concerned  here  with  the  approximation  (5)  in  the  case  which 
includes  the  last  branch  in  (4),  f = 0.  If  f is  a strictly  convex  function^ 
then  it  is  necessarily  non-linear.  However,  the  locus  of  stress  points 
given  by  (5)  as  a function  of  At  is  a straight  line,  so  that  instead  of 
moving  along  the  curved  surface  f = f^,  the  stress  point  moves  along  a 
tangent  to  this  surface.  As  a result,  at  time  t^  + At,  f is  no  longer  equal 
to  f^  but  is  slightly  greater  than  f^. 

This  phenomenon  is  not  just  the  ordinary  one  associated  with  using 
(5)  for  any  nonlinear  time  dependent  problem.  Since  f is  convex,  the 
error  is  always  in  one  direction,  and  if  no  adjustments  are  made,  there 
is  no  limit  to  the  resulting  value  of  f^.  In  the  simple  example  presented 
in  Sec.  5,  an  error  of  up  to  20%  is  easily  produced  in  the  load  for  a 
given  deformation,  and  a solution  is  obtained  for  loads  in  excess  of  an 
upper  bound  on  the  yield-point  load! 

The  above  argument  is  based  on  the  premise  that  f is  strictly  convex. 

However,  if  f is,  instead,  a piecewise  linearyield  function,  then  the  problems 

disappear.  So  long  as  plastic  flow  remains  governed  by  the  same  set  of 

linear  yield  functions,  then  (5)  leads  to  stresses  which  continue  to  satisfy 

f = f^  exactly.  Indeed,  if  we  assume  that  the  boundary  conditions  are  piece- 

wise  linear  in  time,  then  (5)  are  exact  for  a certain  range  of  At,  so  that 

3.  Number  in  brackets  refer  to  the  list  of  references  collected  at  the 
end  of  the  paper. 
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no  error  at  all  is  introduced  by  the  time  integration. 

Now,  the  yield  function  f is  a material  property  which  basically  must 
be  determined  by  experiment.  Experimental  results  will,  of  course,  exhibit 
scatter  due  to  choice  of  sample,  temperature,  previous  material  history, 
etc.,  so  that  the  "true"  f is  not  defined  and  any  f used  in  problem  solving 
is  necessarily  approximate.  Given  this  fact,  it  appears  reasonable  to  choose 
a piecewise  linear  approximation  because  of  its  computational  advantages. 

In  some  simple  problems  where  the  principal  stress  directions  are 
fixed  and  known,  the  piecewise  linear  Tresca  yield  condition  [3]  is  generally 
used  in  preference  to  the  quadratic  Mises  condition  [4],  However,  where 
the  principal  stress  directions  are  unknown  and  vary  with  time,  there  is 
no  generally  accepted  piecewise  linear  yield  condition.  In  this  case,  the 
simplest  procedure  is  to  regard  some  simple  non-linear  criterion  such  as 
the  Mises  yield  condition  as  known  to  within  some  degree  of  accuracy,  and 
then  to  construct  a piecewise  linear  approximation  to  within  the  same 
accuracy. 

In  the  next  two  sections  we  shall  describe,  first  schematically  and  then 
in  formal  detail^ a procedure  for  such  a construction  which  may  be  carried  out 
simultaneously  with  a solution  of  the  boundary  value  problem. 

Sections  4 and  5 will  present  two  simple  examples  to  illustrate  the  method, 
and  the  paper  will  close  with  a discussion  of  the  limitations  of  the  method 
and  some  ideas  for  its  extension. 


2.  Method.  Since  the  proposed  method  is  particularly  useful  in 


L 


obtaining  numerical  solutions  within  the  framework  of  a finite  element 

formulation,  we  shall  present  it  in  those  terms,  but  in  complete  generality 
otherwise.  To  this  end,  consider  a generic  finite  element  k of  the  system 


f 


1 


4. 

whose  stress  configuration  is  described  by  a finite  number  of  generalized 

i< 

efppccpc  Q A 

^a.  Let  the  theoretical  non-linear  yield  condition  for  element  k 
be  where  f^  is  known  only  to  a certain  accuracy.  Then  the  only 

absolute  constraints  on  a theory  are  that  f < and  that  the  element  is 
elastic  for  f < '^q/CI  + “)  where  a defines  the  uncertainty.  For  f between 
these  two  values,  we  are  free  to  construct  any  yield  criterion  which  satisfies 
the  overall  requirement  of  convexity. 

Rather  than  construct  the  detailed  yield  condition  in  advance,  we  shall 
construct  it  in  the  process  of  solving  the  boundary  value  problem.  Before 
giving  details,  we  discuss  the  situation  schematically  in  relation  to  Fig.  1. 

The  problem  is  to  systematically  construct  a convex  polygon  which  will 
lie  in  the  annulus  between  'f'  = 'fp  and  f = fQ/(l  + a).  We  shall  indicate 
four  possible  approaches  and  show  why  one  of  them  is  preferable. 

(1)  initial  tangent  (Fig.  la).  At  the  point  B where  the  stress 
trajectory  first  touches  the  inner  curve,  construct  a tangent  to  the  inner 
curve.  This  method  has  several  drawbacks.  If  the  stress  point  proceeds 
along  the  tangent  to  the  outer  curve  at  point  C,then  some  additional  method 
of  construction  must  be  provided  to  prevent  its  leaving  the  outer  curve.  If 
the  element  unloads  and  reloads  then  a distinction  must  be  made  between  D 
where  we  should  construct  a new  tangent  and  E where  we  would  wish  to  proceed 
to  the  already-constructed  one. 

(b)  external  tangent  (Fig.  lb).  At  the  point  where  the  stress 
trajectory  first  touches  the  outer  curve,  draw  tangents  to  the  inner  curve. 
There  is  a minor  drawback  here  in  that  the  original  trajectory  may  be  outside 


4.  Since  we  will  be  concerned  only  with  the  one  generic  element,  we  will 
suppress  the  superscript  k in  this  section. 


of  the  tangent.  However,  a more  serious  objection  is  that  if  there  are 
more  than  two  generalized  stresses,  then  the  constructed  surface  is  a cone 
which  is  not  piecewise  linear. 

(c)  backtrack  (Fig.  1c).  From  the  point  E where  the  stress  trajectory 
first  touches  the  outer  curve,  backtrack  to  D where  it  most  recently  was 

on  the  inner  curve  and  draw  the  tangent.  This  method  is  some  improvement 
over  the  initial  tangent  one,  but  if  the  loading  now  proceeds  along  the 
tangent  to  F it  is  obvious  that  some  modification  must  be  made. 

(d)  shrink  and  back-track  (fig.  Id).  From  the  point  E with  stress 
coordinates  where  the  trajectory  first  touches  the  outer  curve,  find 
the  point  F with  coordinates  6Q^  on  the  inner  curve  and  draw  the  tangent. 

Then  backtrack  to  G where  the  trajectory  most  recently  crossed  this  newly- 
constructed  tangent.  This  method  works  equally  well  at  H which  was  reached 
along  a tangent. 

Any  back-tracking  method  has  a potential  hazard  in  that  it  is  not  generally 
feasible  to  store  the  entire  trajectory  history.  Suppose,  for  example,  that 
when  the  trajectory  reaches  the  outer  curve  at  W » its  history  is  only 
available  back  as  far  as  M which  violates  the  new  constraint  NP.  In  that 
case,  instead  of  constructing  the  tangent  NP  to  the  inner  curve,  we  construct 
the  parallel  constraint  QR  through  the  point  M. 

On  some  later  loading  the  stress  trajectory  may  encounter  a previously 
drawn  tangent,  in  which  case  plastic  flow  will  be  governed  by  it.  For 
example,  along  TU  the  trajectory  cannot  continue  to  the  outer  curve  to  trigger 
a new  tangent.  Rather,  it  will  move  along  the  already  established  tangent 
HU. 

It  appears  from  this  schematic  analysis  that  the  shrink  and  back-track 
method  does  not  suffer  from  the  drawbacks  of  the  other  three,  so  we  shall  now 


r 


proceed  to  formalize  the  procedure. 


3.  Procedure.  We  assume  without  any  real  loss  in  generality  that  the 
given  loads  and  any  prescribed  displacements  are  all  piecewise  linear  in  time. 
If  they  are  first  given  by  explicit  non-linear  equations  these  equations 
can  be  replaced  by  piecewise  linear  approximations  before  starting  the 
problem.  Also,  we  take  the  finite  element  model  as  given. 

b 

Let  denote  the  ath  generalized  stress  in  element  k,  and  let  Q be 

l( 

the  array  of  all  Q . Similarly,  let  q , p,  and  u denote,  respectively, the 
arv'-'vs  of  generalized  total  strain,  generalized  plastic  strain,  and  generalized 
acement. 

Ihe  total  and  plastic  strain  rates  are  related  by 

q = £9  + p (6) 

where  c is  the  array  of  elastic  constants.  If  stress,  strain,  and  plastic 
strain  are  all  initially  zero,  (6)  may  be  integrated  to  yield 

g = cQ  + p (7) 

and,  for  the  solution  of  any  problem,  u can  be  determined  from  q . Therefore, 
the  rate  problem  need  be  concerned  only  with  the  stress  rates  Q and  plastic 
strain  rates  p. 

For  element  k,  plastic  flow  will  be  zero  or  will  be  governed  by 
one  or  more  linear  yield  flats 


J aj  ^cx 


(8) 


where  the  summation  convention  applies  only  to  Greek  subscripts.  At  any 
stage  in  the  solution  procedure,  a finite  number  of  the  L.  will  be  defined, 

J 

and  at  the  start  of  the  stage  each  L-  will  either  be  equal  to  one  or  less 

J 

k k 

than  one.  Flats  with  L.  < 1 will  be  inactive  and  flats  with  L.  may  be  either 

J J 

active  with 

l'^  s a*^.  (5*^  = 0 (9) 

J aj  a 


(9) 


7. 


or  inactive  if 


= a^.  < 0 

J aj  a 


The  plastic  strain  rates  in  element  k are  given  by 


= z:.  a*",  x'l 


AT  > 0 


(10) 


(11) 


J J aj  j ’ 

where  the  summation  is  to  be  taken  only  over  the  active  flats  in  element  k. 

Equation  (6)  may  be  substituted  in  the  compatibility  equations  to 
express  them  as  linear  relations  between  Q and  p,  and  the  p may  be  replaced 
by  the  A from  (11).  Each  active  flat  will  provide  an  equation  in  the  form 
(9)  and  introduce  an  unknown  plastic  multiplier  Therefore,  for  a properly 
formulated  boundary  value  problem,  equilibrium,  compatibility, and  (9)  will 
provide  sufficient  linear  equations  to  determine  all  unknown  rates.  Although 
the  question  of  determining  which  flats  are  active  may  be  a vexing  one,  we 
shall  assume  that  it  can  be  satisfactorily  resolved  and  regard  the  rate 
problem  as  solvable. 

The  total  problem  will  consist  of  a finite  number  of  stages.  In  each 
stage,  the  rate  problem  will  first  be  solved  using  given  values  at  the 
beginning  of  the  stage.  Next,  a number  At  will  be  obtained  which  denotes 
the  duration  of  the  stage,  and  possibly,  a new  flat  may  be  created  for 
future  stages.  Finally,  values  at  the  end  of  the  stage  will  be  computed  to 
act  as  initial  values  for  the  next  stage. 

Figure  2 gives  a gross  outline  of  the  process.  We  begin  with  ini- 
tial values  for  Q,  p,  and  t.  Normally  these  will  be  zero,  but  it 
might  be  desirable  to  begin  in  the  middle  of  a previously  solved 
probl em . 

The  rate  problem  is  solved  including  verification  that  all 

k 

non-zero  A*^  are  positive  and  that  all  inactive  flats  where  = 1 
J 

satisfy  (10) , 

Next,  three  numbers  are  determined.  At^  is  the  time  remaining  unti 


8. 


some  exterior  requirement  forces  an  end  to  the  stage.  This  will  normally 

be  the  end  of  a linear  interval  in  the  initial  piecewise  linearization  of 

applied  loads  and  displacements,  but  it  may  be  presented  in  other  terms. 

Atp  is  defined  as  the  minimum  positive  time  interval  over  all  elements  k for 

k k°  xk 

the  element  to  reach  its  nonlinear  yield  value:  f (Q^  + At  Q^)  = 1 
corresponding  to  the  outer  curve  in  Fig.  1.  Finally,  At^  is  defined  as 
the  minimum  positive  time  interval  over  all  elements  k and  over  all  inactive 
flats  for  the  element  to  reach  the  limiting  flat  value  + AtO^)  = 1. 

J J cx  ot 

Let  the  smallest  of  these  three  intervals  be  denoted  by  At*.  If  At*  is 

either  At^  or  At^,  i.e.,  if  the  stage  is  terminated  by  either  an  external 

k k 

cause  or  a previously  defined  inactive  flat  so  that  all  f (Q^)  < 1,  then 
we  set  At  = At*  and  proceed  immediately  to  ending  the  stage.  However,  if 
At*  = Atp,  then  it  is  necessary  to  form  a new  flat  and  back-track.  To  this 
end,  let 

+ At*  Q (12) 

a a a ' ^ 

be  the  stress  point  responsible  for  terminating  the  interval  (e.g.,  point 
E in  Fig.  Id),  and  define  the  stress  point  by 

Q^'  = m*  f(Q^’)  = l/(l+a)  (13) 

When  f is  homogeneous  in  the  generalized  stresses,  it  may  be  more 

k k 

meaningful  to  define  the  inner  yield  surface  by  f [0-e)Q^]  = 1 where  e 

represents  the  estimated  error  range  in  the  yield  stress.  If  f is  homogeneous 

of  degree  n,  then  the  two  parameters  are  related  by 

e=l-(l+a)*^  (14) 

Next  a new  flat  is  constructed  for  element  k tangent  to  the  inner 
k ' 

surface  at  Q : 

a 

(q'^-q'^')  — . (Q*^')  = 0 


(15) 


9. 


Defining 


k ^ a 

Qg3f/aQg 


we  can  rewrite  (15)  in  the  form  (8).  Here  j is  to  be  assigned  the  value 

I 

j'  + 1 where  j'  was  the  number  of  the  last  flat  to  be  constructed  for  element 


The  interval  At  is  now  defined  by 

so  that  the  stress  point  lies  on  the  new  flat  fe.g.  point  G on  trajectory 
DE  in  Fig.  Id).  If  At  is  positive  it  denotes  the  end  o'*"  the  stage  and  we 
compute  the  final  values.  However,  if  At  is  negative  would  have  to  back- 
track into  the  previous  stage,  a process  which  is  not  generally  desirable. 

In  this  case,  we  cancel  out  the  current  stage  by  setting  At  = 0,  but  we  do 
add  a new  flat  which  will  necessarily  be  active  and  hence  change  the  rate 

solution.  The  new  flat  is  to  be  parallel  to  (15)  which  led  to  the  negative 

k° 

At,  and  is  to  pass  through  the  point  where  the  stage  initiates.  Thus, 

in  place  of  (16)  we  define  i, 

. Of/3Q  )'  , , 

» r “-r- 

Qg  (Sf/3Qg)' 

where  the  gradients  are  still  to  be  evaluated  at  the  original  point  of  inner 
tangency  defined  by  (13). 

With  At  determined  under  all  contingencies,  the  stage  is  ended  with 

t = t^  ^ At 

Q = Q°  + At  Q 


p = p + At  p„ 

a a 


10. 


strains  are  computed  from  (7),  and  displacements  from  the  strain-displacement 
relations.  Any  desired  output  is  generated  and  the  new  state  is  tested 
to  see  if  the  problem  is  completed.  If  not,  the  left-hand  sides  of  (19) 
are  taken  as  the  initial  values  for  the  next  stage,  and  the  process  is 
repeated. 

4.  Frame  example.  As  a first  example  we  consider  the  frame  shown  in 
Fig.  3a.  The  top  and  bottom  bars  are  rigid, and  the  top  bar  is  fixed  so 
that  the  deformation  has  three  degrees  of  freedom  as  shown  in  Fig.  3c; 

Fig.  3b  gives  the  sign  convention  for  the  stress  resultants. 

We  consider  a simplified  model  in  which  the  deformation  is  confined 
to  hinges  at  the  ends  of  the  bars.  Then  so  long  as  at  least  one  bar 
remains  elastic,  which  turns  out  to  be  the  case  up  until  collapse,  the 
shear  forces  S^.  are  all  zero,  the  moments  and  normal  forces  at  either  end 
of  each  bar  are  equal,  and  Ln  = cf/Z  so  that  the  frame  has  only  two  degrees 
of  freedom.  Further,  since  conditions  in  the  hinges  at  either  end  of  a beam 
are  the  same,  we  can  carry  out  the  analysis  in  terms  of  three  elements  each 
representing  the  hinges  in  one  beam. 

Let  and  be  the  yield  moment  and  force, 'respectively,  and  define 
dimensionless  generalized  stresses  and  a dimensionless  load  by 

m.  = M./M  , n.  = N./N  , p = P/N  (20) 

1 r 0 1 10^  0 

Then  moment  and  vertical  equilibrium  of  the  lower  bar  lead  to 

3n.|  + n^-n^  = h(m.|+m2+m2)  (21a) 

n^+n^+n^  = P (21b) 

where  h = 

0 0 

Generalized  strains  rorrespondi ng  to  m^.  and  n^  are  the  slope  discontinuity 
0^.  across  the  hinge  and 
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" 0/h)(u/L) 


(22) 


where  u^.  is  the  axial  length  discontinuity.  In  terms  of  the  kinematic 


variables  4>  and  6,  the  generalized  strains  are  given  by 


2he-j  = 6-3(J) 
2he  = 6+(j) 


2he2  = 

01  = 02  = 03  = -(p/2 


(23) 


In  the  numerical  example  to  follow  we  shall  take  h = 0.1. 

As  described  in  [5]  the  elastic  constitutive  equations  for  this  simple 
frame  model  can  only  be  expressed  in  terms  of  an  indeterminate  magnitude 
constant  which  represents  a scale  factor  for  the  displacement  configuration. 
Without  loss  of  generality,  we  take  this  constant  to  be  unity,  whence 
with  h = 0.1  the  elastic  rate  equations  can  be  written 

n-j  = 2.5(6-341)  n2  = 2.5(6-41) 

n^  = 2.5(6+^)  m^  = m2  = m^  = -4/3 


(24) 


In  order  to  illustrate  a variety  of  loading  behavior,  we  shall  choose 
the  yield  stresses  shown  in  circles  in  Fig.  3a  for  the  three  bars  as  in 

Drucker's  truss  [6].  Then  the  yield  curves  each  consist  of  two  parabolas 

.2 

I III/  ZU  = I 

(25) 


f^  = n^  ± m = 1 f2  = (n/20)^  ± m/20  = 1 


f^  = (n/12)^  ± m/12  = 1 

Finally,  we  shall  choose  the  factor  a between  the  inner  and  outer  curves  to 
be  0.1. 

If  element  i is  on  a single  yield  flat 


L^.  = an  j + bm^  = 1 


(26) 


then  the  constitutive  equations  can  be  explicitly  solved  for  the  generalized 
stress  rates  to  yield 
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ii.  = b(2be^.+a4i)/(3a^+4b^)  = -an^/b  (27) 

and  thR  condition  for  continuing  plastic  flow  derived  from  X > 0 is 

3ae.-2bi  > 0 (28) 

When  element  i stays  on  two  independent  flats  of  the  form  (26), 
the  only  solution  is  n^.  = m^  = 0 in  place  of  (27).  In  this  case,  the 
inequality  (28)  must  be  satisfied  for  both  flats. 

We  consider  the  behavior  of  the  frame  as  P is  slowly  increased  until 
collapse  occurs.  It  is  convenient  to  take  6 (which  is  monotonical ly  increasing 
under  increasing  P)  as  the  time  variable,  whence  the  load  is  determined  by 
Eq.  (21b). 

All  elements  are  elastic  in  the  first  stage,  which  ends  at  6^  = 1.963 
when  f^  = 1 with  (n^ , m^ ) = (0.907,  -0.178).  The  corresponding  shrunk  point 
is  B(n^,  m.|)  where  B is  the  solution  of 

(0.9076)^  + (0.178B)  = 1/1.1  (29) 

which  leads  to  (n^ , mp  = (0.860,  -0.169).  The  tangent  line  to  the  inner 
yield  curve  at  this  point  is 

= 1 .043n^-0.606m^  = 1 (30a) 

and,  repeating  the  elastic  solution,  we  see  that  the  limiting  value  is  6 = 1.863. 
Final  elastic  values  for  stage  1,  which  will  serve  as  initial  values  for  stage 
2,  are  recorded  on  line  1 of  Table  1. 

In  stage  2,  element  1 proceeds  along  until  it  hits  the  other  yield 
curve  f|  = 1 at  A6p  = 0.394  with  (n.j  , m.j ) = (0.934,  0.036).  Shrinking  then 
leads  to  the  new  tangent  line 

= 1.049n.j  + 0.561m^  = 1 (30b) 

with  A6  = 0.313  and  the  final  values  recorded  in  line  2 of  Table  2.  It 
is  interesting  to  note  that  the  end  value  of  m^  is  negacive  even  though 
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= 1 was  the  limiting  yield  curve. 

During  stage  3 element  1 remains  at  the  intersection  of  L-|.|  and  L-|2. 
and  the  stage  limit  is  determined  by  element  3 reaching  f ^ = 1 . This  leads 
to  the  tangent 

13^  = 0.0874n3-0.0462m3  = 1 (30c) 

and  the  values  in  line  3 of  Table  1. 

In  stage  4 the  frame  exhibits  behavior  similar  to  Drucker's  truss  [6] 
in  that  element  1 leaves  both  linear  yield  faces  and  unloads  elastically. 

A6p  is  determined  by  f^  = 1,  and 

L32  = 0.0874n3  + 0. 0469m 3 = 1 (30d) 

For  all  further  loading,  element  3 remains  at  the  point  (n3,  m3)  = 
(11.443,  0.003)  which  is  the  intersection  of  L3.1  and  L32.  At  the  end  of 
stage  5,  during  which  it  is  elastic,  element  1 yields  in  compression  but 
with  a substantial  negative  moment:  (n^,  m^  = (-.677,  -.452).  During  stage 
6 it  yields  on 

L^3  = -0.995n^  - 0.724m^  = 1 (30e) 

and  during  stage  7 it  leaves  L-|3  and  moves  along 

= -1.049n^  - 0.557m^  =1  (30f) 

Stage  7 terminates  on  the  tangent 

= -1.057n-|  + 0.563m^  = 1 (30g) 

based  on  f|  = 1,  and  from  here  on  element  1 remains  at  the  intersection  of 

*■14  *‘15' 

With  both  bars  1 and  3 capable  of  two  degrees  of  freedom  in  their  strain 
rates,  the  frame  can  carry  additional  load  only  because  of  bending  strength 
in  bar  2.  After  a very  small  increment  of  load  accompanied  by  large  increases 
of  6 and  4>,  bar  2 is  restricted  by  f2  = 1 and  the  final  tangent  plane 

1-2^  = 0.0503n2-0.0353m2  = 1 (30h) 
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If  the  load  is  maintained  at  the  value  p = 25.484,  then  a solution  exists 
with  all  stress  rates  equal  to  zero  and  6 = 4>  being  any  positive  number. 
Therefore,  collapse  occurs  with  element  1 on  L-j^  and  element  2 at 
a fixed  point  on  L21,  and  element  3 on  L^-]  and  122-  Final  collapse  values 
are  given  on  line  8 of  Table  1.  Figure  4 shows  the  stress  trajectories  and 
yield  lines  of  the  three  elements,  and  Fig.  5 shows  the  load-deformation 
diagrams. 


5.  Plane  stress  example.  To  show  that  the  method  works  equally 
simply  when  there  are  more  than  two  generalized  stresses  per  element,  we 
consider  the  simple  plane  stress  example  shown  in  Fig.  6.  The  triangular 
plate  is  fixed  along  its  vertical  edge  AB  and  an  arbitrary  in-plane  load 
is  applied  at  C.  We  divide  the  plate  into  four  triangles  as  shown  and 
assume  a linear  displacement  field  in  each  so  that  there  are  six  degrees  of 
kinematic  freedom  described  by  the  nodal  displacements.  In  order  to  carry 
out  the  subsequent  analysis  in  dimensionless  terms,  we  shall  take  these  to 


(u.,v^=  (2E/La^)(U.,V.)  i = 1,  2,  3 (31 

where  U^,  are  the  physical  displacements. 

4 3 

Evidently  the  strains  in  triangles  3 and  4 are  all  equal  and  = 0, 

so  that  there  will  be  eight  generalized  strai  ns,  def  i ned  by 

2qi  = 2EeX  = = V2-V3 

2^3  ^ 2Ey‘/o^  = U2-U3  + v^-V3 


2q.  = 2Eet/a„  = u- 


2^5  = 2^^y%  = "2-  "3 


2^6  = 2EY,y/a^  = U2-U3+V3 


P.Eeyo, 
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The  corresponding  generalized  stresses  are  the  average  dimensionless  stresses 
in  each  triangle.  Since  constant  strains  lead  to  constant  stresses  for  an 
elastic/perfectly-plastic  material,  these  are  simply 


Q,  = a,,  t^.  cl,  ^3^) 


(33) 


Although  is  not  necessarily  zero,  it  does  no  work  and  hence  is  not  a 

generalized  stress;  it  can  always  be  eliminated  from  the  problem  by  the 

3 ... 

constraint  = 0 on  the  constitutive  equations. 

The  eight  Eqs.  (32)  may  be  inverted  to  yield  explicit  expressions  for 

displacements  in  terms  of  strain,  and  two  compatibility  relations: 


^2  ^7  '^2  ^^3  ^*^4  ''3  (34a) 

= ^(q^+q^)  V.,  = 2(q3+2q^+qg-q7)  (34b) 

q2-qg  = 0 2(q^+q5+qg)-(q7+qg)  =0  (34c) 

3 

For  a state  of  plane  stress  with  v = 1/3  and  = 0,  the  constitutive 
equations  (7)  may  be  written 

q^  = Q^-Q2/3+p^  q2  = Q2-Q^/3+P2  q3  = 8Q3/3+P3 

q^  ~ Q4"Qg/3'''P4  q^  ~ ^5”^4^^^P5  (35) 

q^  = 16Qy/9+Py  qg  = 16Qg/3+Pg 

Finally,  the  equilibrium  relations  between  the  generalized  stresses  are 


Q3+Qg+2Q7  = 0 Q2+Qg+2Qg  = 0 

Q,+Q3-Q4*Q5  ■ 0 ■ 0 


(36) 
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= = (37) 

We  consider  a problem  in  which  the  tip  displacement  u.j  = u,  v-j  = v 

is  prescribed  as  shown  by  the  heavy  curve  1-7  in  Fig.  7.  Point  7 is  determined 

by  the  condition  P = 0,  and  at  this  point  P will  be  negative.  The  problem 

y ^ 

is  concluded  by  increasing  P^^,  through  zero,  until  collapse  occurs. 

Since,  for  most  of  the  history,  u and  v are  prescribed,  it  is  convenient 
to  use  these  as  the  independent  time-like  variables.  To  this  end,  Eqs. 

C35)  are  first  substituted  in  (34)  and  the  results  are  combined  with 
(36)  to  provide  eight  equations  for  the  generalized  stresses  in  terms  of  u, 

V,  and  the  plastic  strains.  It  turns  out  that  triangle  3 is  always  elastic, 
so  that  P^  = Pg  = 0 at  all  times,  The  resulting  equations,  written  in  terms 
of  rates,  are  displayed  in  Table  2. 

The  Mises  yield  condition  is  taken  for  the  outer  yield  surface,  hence 

fl  = Qi-QiQ2+Q^3Q^ 

^2  " (38) 

fj  = (7/9)Q^  + 3q2 

3 

In  writing  f.^  we  have  used  the  elastic  relation  c = a - Q«/3  = 0 to 

y y ^ 

eliminate  a^,  since  triangle  3 never  becomes  plastic.  The  inner  curve  will 
be  defined  by  a yield  stress  95%  of  o^,  hence  it  follows  from  (14)  that 
a = .02598. 

We  identify  eight  stages  by  the  numbers  in  Fig.  7 and  use  decimal  sub- 
stages to  indicate  different  plastic  behavior.  In  stage  1,  u = 0,  v is 
decreased  to  - 2.25,  and  all  elements  are  elastic.  Therefore,  all  p^.  = 0 
and  Table  2 immediately  gives  the  stresses  in  terms  of  v.  The  stage 
terminates  externally  when  -v  = 2.25. 
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Stage  2 is  initially  elastic  as  is  increased  with  held  constant, 
but  at  Up=  2.221,  = 1 . with  a*  = (Q^  Q2.  Q3)  = (0.847,  0.204,  -0.372). 

The  associated  stress  point  on  the  inner  curve  is  o'  = 0.95o*  and  the  tangent 
plane  through  o'  is 

= 0.784Q^-0.231Q2-1.174Q3  = 1 (39) 

Backtracking  the  elastic  solution,  we  find  that  stage  2.1  ends  with 
u=2.055  with  the  stress  values  given  on  line  2.1  of  Table  3. 

In  stage  2.2,  triangle  1 is  on  L-ji,  hence  it  follows  from  (11)  and 
(39)  that 


V3^,P2,P3)  = A(0. 784, -0.231  , -1.174) 


(40) 


Substitution  of  (40)  in  Table  2 and  the  results  in  (8)  for  L-ji  leads  to 

A = 0.375  u - 0.182V  (41) 

Equations  (40)  and  (41)  with  v = 0 can  then  be  substituted  in  Table  2 to 
give  explicit  expressions  for  the  stress  rates.  This  stage  terminates  when 
the  displacement  prescription  changes  at  AUp  = 0.195,  with  results  given  on 
the  lines  2.2  in  Tables  3 and  4. 

Stage  3 with  u * v is  also  divided  into  two  parts.  In  stage  3.1, 
triangle  1 remains  on  L-ji  until  f^  = 1 again  for  AUp  = 1.511  which  leads  to 
the  new  tangent  plane 

L^2  = 0.9640^-0. 36702-0. 697Q3  = 1 (42a) 

and  backtracking  shows  that  the  stage  ends  at  Au  = 0.730.  In  stage  3.2, 
triangle  1 leaves  but  stays  on  1^2  until  stage  3 terminates  with 
AUp  = 1.520. 

Stage  4,  in  which  v is  increased  to  4.5  with  u held  constant  at  4.5, 
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begins  with  triangle  1 leaving  L-|2.  so  that  the  entire  plate  is  elastic. 

Now  f2  = 1 becomes  the  violated  condition,  and  in  stage  4.2  triangle  2 
is  on  the  tangent  plane 

L21  = 1 .030Q^-0.679Q^-0.185Qg  = 1 (42b) 

In  stage  5,  v increases  from  4.5  to  9.0  and  u = -v.  There  are  four 
substages  as  triangle  2 unloads  elastically,  but  triangle  1 repeatedly 
encounters  f-j  = 1.  Thus,  stage  5.1  is  elastic,  and  in  stages  5.2,  5.3, 
and  5.4,  triangle  1 is  respectively,  on 


^3  = 0.114Q^-0.287Q2+1.753Q3  = 1 

(42c) 

^4  = -0.211Q^-0.065Q2+1.753Q3  = 1 

(42d) 

^3  = -0.521Q.|+0.136Q2+1.564Q3  = 1 

(42e) 

Thus  far  we  have  not  mentioned  symmetry  or  isotropy.  Whenever  one 
yield  flat  is  formed,  one  may  wish  to  immediately  introduce  various  symmetric 
linear  constraints.  Thus,  if  the  loading  has  established 

L = aa^+bOy+CT  =1  (43) 

as  a linear  yield  condition,  one  has  the  option  of  immediately  defining 
any  of 


L 

= -aa^-bOy-CT  = 1 

(44a) 

L"  = aa +bo  -CT  =1 
X y 

(44b) 

L"’ 

= -ao^+bOy+CT  = 1 

(44c) 

L " " 

= bo  +ao  +CT  = 1 
X y 

(44d) 

representing,  respectively,  complete  symmetry  in  tension-compression,  shear 
symmetry,  symmetry  of  in  tension-compression,  and  isotropy.  If  more  than 


one  of  (44)  are  used  they  will  obviously  combine  to  give  further  combinations. 

To  illustrate  the  application  of  this  idea,  we  shall  adopt  only  (44a) 
so  that  for  each  yield  flat  defined  by  (30)  or  (42),  another  flat 
Lj^j  =-  L|^j  is  also  defined. 

These  additional  yield  flats  first  become  relevant  during  stage  6 in 
which  V = 9.0  and  t = -u  is  increased  from  0 to  9.0.  In  stage  6.1,  triangle 
1 remains  on  L.|g  and  at  Atp  = 2.275,  f-j  = 1 again.  However,  at  At|^  = 1.336, 
triangle  1 encounters  flat  i.e.,  L-j-j  as  given  by  Eq.  (39)  takes  on  the 
value  -1.  Stages  6.2  and  6.3  find  triangle  1 on  L^.j  and  respectively, 
but  then  the  surface  f^  = 1 is  encountered  again  and  a new  flat 


= -0.8870^+0. O53Q2+O.595Q3  = 1 


is  defined. 


In  stage  6.4,  triangle  1 remains  on  both  1^3  and  L-jg,  but  since  three 
stresses  are  involved  this  does  not  necessarily  correspond  to  a constant 
stress.  Instead,  the  plastic  strain  rates  must  be  written 

(-0.964,0.367,0.697) 

(46) 

+ A2(-0. 887, 0.053, 0.595) 

and  the  two  additional  equations  L.jg  = 1^2  ^ 0 solved  simultaneously  to  yield 

A^  = 0.033t  A2  = 0.421t  (47) 


Substitution  of  (46)  and  (47)  in  Table  2 then  leads  to  the  stress  rates. 

Stage  6.4  ends  when  triangle  2 encounters  flat  L^-j-  In  the  next  stage, 
Llo.  Lie,  and  LI,,  are  all  active.  In  a load-controlled  test,  activation  of 

\ C I 0 c\ 

three  yield  planes  would  generally  cause  collapse  with  all  0.j  = 0>  but  this 
solution  is  not  consistent  with  v = 0.  Instead,  (46)  must  be  used,  together 


20. 


wi  th 


(p^,Pg,Pg)  = ^(-l. 030, 0.679, 0.185)  (48) 

and  1^2  “ *-i6  ^ '*"21  ~ ^ obtain 

= 0.256t  A2  = 0.122t  p = 0.143t  (49) 

whence  Table  2 again  provides  the  stress  rates.  Stage  6.5  terminates  when 
u = -9.0. 


In  stage  7,  u = v = -t.  It  turns  out  that  under  this  loading  triangle 
1 leaves  L^g  but  remains  on  1^2’  triangle  2 remains  on  L^-j . This  stage 
is  to  continue  until  = Qj/2  = 0 which  leads  to  At^  = 3.926.  However, 
triangle  1 hits  a new  portion  of  f.|  = 1 at  Atp  = 2.692,  whence  stage  7.1 
ends  with  formation  of 


=-  1.047Q^  + O.47OQ2  + O.I53Q2  = 1 (50) 

and  At  = 0.937.  Flats  L.|y  and  L^-]  are  the  active  ones  in  stage  7.2  and 

this  time  the  external  condition  Q3  = 0 governs  the  end  of  the  stage. 

Stage  8 is  a load  controlled  one,  so  it  is  convenient  to  solve  the  first 

and  third  equations  in  Table  2 for  u and  v in  terms  of  P (and  P = 0). 

A y 

When  the  results  are  substituted  in  the  remaining  equations  in  Table  2,  explicit 
dependence  cn  p^  and  p^  disappears,  and  the  resulting  stress  equations  can 
be  simply  written  in  terms  of  the  parameters 

t = 0.08ya^ 

r = -0.06( P4+Pg+Pg)+0.52(pg-p2)  (51) 

s = -0.18(p4-tp*5'*'P5)+0-0^(P5-P2^ 


vHPiMBiapppn 
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the  results  being 

Q-|  = 25t  ^2  ^ 

= 22t+s  O5  = 2t-r+s  Qg  = -3t+s  (52) 

O7  - Os  » -0/2)Q6 

All  triangles  are  elastic  in  stage  8.1,  which  terminates  when  triangle 
1 encounters  and  remains  there  until  collapse.  Triangle  2 is  still 
elastic  in  stage  8.2,  but  returns  to  L21 , (Eq.  42b)  for  stage  8.3.  The  limiting 
condition  for  stage  8.3  is  f2  = 1 which  leads  to  the  final  plane 

L22  = 0.932Q^-0.684Qg-0.725Qg  = 1 (53) 

Collapse  now  occurs,  since  Q^.  = 0 is  the  unique  solution  of  Eqs.  (52)  and 
~ ^21  ^ ^22  ~ three  plastic  multipliers  are  positive. 

At  each  substage  of  the  solution, the  plastic  strain  rates  are  ob- 
tained, and  a running  record  is  kept  of  their  accumulated  values  in  Ta- 
Dle  4.  However,  Eqs.  (34)  and  (35)  give  the  displacements  and  strains  in 
finite  form, so  that  they  are  easily  computed  from  Tables  3 and  4 whenever 
they  are  desired.  During  the  def  ormat  i on -con  t ro  1 1 ed  process,  the  loads  are 
given  by  (37),  and  when  loads  are  controlled  t ip  def  0 rma  t i on  s are  given  by  (34b  ). 

It  is  of  some  interest  to  compare  the  preceding  solution  with  an  unadjusted 
use  of  the  Mises  yield  condition  using  approximately  the  same  number  of  time 
steps.  In  this  analysis,  the  limit  and  end  of  a step  is  either  At^,or  Atp 
when  an  elastic  triangle  reaches  the  current  yield  value  f^.  When  f?  has 
been  reached  plastic  flow  will  take  place  along  the  tangent  plane  to  the 
surface  f^  = f?  so  that  at  the  end  of  the  step  the  value  of  f?  must  be  recomputed 
as  the  then  current  value  of  f . For  a fair  comparison,  stages  2 and  4 each 
consist  of  one  elastic  and  one  plastic  step,  stage  3 is  divided  into  two  steps 
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with  Au  = 1.125  and  stage  5 consists  of  an  elastic  step  followed  by  three 
plastic  steps  with  equal  Av.  Stage  6 is  first  subdivided  into  five  steps 
with  -Au  = 1.80  in  each,  but  the  last  step  is  further  split  in  two  by  triangle 
2 becoming  plastic.  Stage  7 is  taken  in  equal  steps  of  Au  = -2.5,  and  it  is 
found  that  Py  = 0 part  way  through  the  second  step.  The  first  step  in  stage 
8 is  elastic  until  triangle  2 becomes  plastic.  Step  8.2  is  terminated  by 
triangle  1 becoming  plastic  and  subsequent  steps  are  chosen  so  that  AA,  the 
integrated  plastic  strain  multiplier  for  triangle  1 is  equal  to  1.00. 

It  is  found  that  the  load  does  not  increase  significantly  after  the  end  of 
8.2,  but  a collapse  mechanism  occurs  only  at  the  end  of  stage  8.5.  Values 
of  stress  and  plastic  strain  according  to  this  analysis  are  presented  in 
Tables  5 and  6.  Table  6 also  shows  the  current  yield  values.  By  the  end  of 
the  loading,  because  of  finite  motion  on  tangent  planes,  the  yield  value  for 
triangle  2 has  increased  10%,  and  for  triangle  it  has  increased  67%. 

Figure  7 gives  a pictorial  history  of  the  load  and  tip  deformation. 

The  deformation  histories  are  prescribed  through  the  first  7 stages,  but 
differ  noticeably  during  the  final  increase  of  P^,  wheras  the  load  histories 
gradually  diverge  through  the  first  7 stages  and  show  a 25%  difference  in 
the  predicted  collapse  loads.  The  light  ellipse  in  Fig.  7 is  an  upper  bound 
on  the  yield-point  interaction  curve,  corresponding  to  plastic  deformation 
of  only  triangle  1.  The  unadjusted  analysis  predicts  loads  in  excess  of 
collapse  from  stage  5.1  through  stage  7,  and  also  for  the  final  collapse 
load. 

Figure  8 presents  the  same  information  from  a different  viewpoint  as  a 
load-deformation  diagram  in  the  horizontal  direction.  Figure  9 gives  an 
exaggerated  view  of  the  total  deformation  at  the  end  of  stage  6.  Although 
the  tip  deformation,  being  prescribed,  is  the  same  for  both  theories,  the 
strains  and  other  nodal  deformations  are  measurably  different  between  the  two. 
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6.  Conci usions  A method  has  been  described  and  illustrated  for  automatic 
piecewise  linearization  of  finite  element  analyses  for  1 inearly-elastic/ 
perfectly-plastic  structures.  According  to  the  method,  material  behavior 
is  originally  described  by  a nonlinear  outer  yield  surface  and  a related  inner 
yield  surface.  As  the  loading  progresses,  yield  flats  are  formed  which  are 
linear  constraints  lying  wholly  between  the  two  surfaces.  Therefore,  at 
any  instant,  the  current  yield  condition  consists  of  any  number  of  yield 
flats  plus,  possibly,  a portion  of  the  original  outer  surface.  Although 
the  current  condition  may  change,  at  any  instant  it  is  fully  contained  within 
the  yield  condition  at  any  previous  instant. 

The  examples  were  purposely  chosen  to  be  very  small  in  order  to  clearly 
illustrate  the  principles  involved.  In  implementing  the  procedure  for 
large  scale  problems,  various  refinements  might  be  considered. 

The  use  of  symmetry  and  isotropy  to  more  rapidly  create  linear  yield 
flats  has  already  been  commented  on  in  Section  5.  It  should  be  pointed 
out  that  isotropy  is  used  in  a very  limited  sense.  If  the  principal  stress 
directions  are  variable,  then  no  piecewise  linear  yield  function  can  be 
fully  isotropic,  since  a general  rotation  of  axes  transforms  a linear 
function  into  a non-linear  one.  The  isotropy  referred  to  amounts  only  to 
interchange  of  fixed  coordinate  axes  or,  equivalently,  to  rotations  which 
are  multiples  of  90°. 

Another  possible  refinement  is  to  use  the  same  piecewise  linearization 
with  respect  to  all  elements.  This  would  obviously  reduce  storage  requirements. 
A minor  problem  would  exist  if  the  stress  point  for  element  2,  say,  was 
between  the  inner  and  outer  yield  curves,  but  fell  outside  of  a yield  flat 
newly  created  by  element  1.  However,  this  aspect  could  easily  be  handled 
by  a method  similar  to  the  loss  of  ancient  history  problem  as  at  point 
M in  Fig.  Id. 
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In  order  to  reduce  the  number  of  time  steps,  it  might  be  desirable  to 
examine  all  elements  v^ien  a stage  is  terminated,  and  to  construct  yield 
flats  for  any  whose  stress  points  are  outside  the  inner  curve,  even  if  they 
are  not  yet  on  the  outer  curve. 

Two  drawbacks  to  the  procedure  should  be  noted.  First,  if  two  or  more 
problems  are  solved  for  the  same  structure  as  part  of  the  same  program  so 
that  each  problem  starts  with  any  yield  flats  developed  from  previous  ones, 
then  the  solutions  may  differ  depending  upon  the  order  of  solution.  Second, 
a stress  point  strictly  between  the  two  curves  may  be  reached  early  in  the 
loading  history  which  is  outside  of  a flat  established  later.  It  is  even 
possible  that  such  a point  could  correspond  to  plastic  collapse  of  the 
structure  according  to  the  final  yield  criterion.  Although  either  of  these 
situations  is  theoretically  embarrassing  they  should  not  prove  to  be  of 
great  practical  significance,  since  the  current  yield  condition  is  always 
bounded  by  the  original  inner  and  outer  surfaces  which  can  be  chosen  to  be 
reasonably  close  to  each  other. 

The  method  is  readily  adaptable  to  strain-hardening,  and,  indeed,  might 
be  somewhat  simpler  to  apply.  For  a strain-hardening  material,  let  the 
initial  yield  condition  for  an  elastic  element  be  the  inner  yield  curve,  and 
terminate  the  elastic  stage  when  this  limit  is  reached  as  at  point  C in 
Fig.  10. 

Construct  the  yield  flat  JCD  and  assume  loading  continues  on  this  flat 
according  to  whatever  hardening  law  is  being  used.  During  this  stage  the 
current  yield  condition  consists  of  the  flat  JCD  together  with  that  portion 
of  the  outer  yield  curve  which  has  not  been  cut  off  by  JCD.  If,  for  example, 
the  hardening  is  isotropic,  this  entire  configuration  would  expand  as  the 
stress  point  moved  outside  the  original  flat  JCD.  This  phenomenon  is  easily 
incorporated  into  Fig.  10  by  letting  the  stress  scale  shrink  to  account  for 


/ 


¥ 


I 

! 

1 

1 

i 


! 

I 


25. 


the  hardening.  Similarly,  if  the  hardening  is  kinematic,  the  stress  origin 
can  be  moved. 

If  the  stress  point  reaches  D (in  its  new,  hardened  position),  then  the 
shrink  and  backtrack  technique  is  used  just  as  in  perfectly  plasticity  to 
form  a new  flat  ENK. 

If  unloading  occurs  at  point  F,  the  behavior  will  be  elastic  and  the 
current  yield  condition  consists  of  the  existing  flats  JCD  and  ENK,  the 
inner  yield  curve  QIT  which  is  not  influenced  by  the  flats,  and  conical 
portions  such  as  AQK  connecting  these  two.  Note  that  this  surface  is  not 
convex. 

Plastic  flow  may  recommence  in  three  different  manners.  For  trajectory 
FGI  it  will  start  with  formation  of  a new  tangent  flat  at  I.  For  trajectory 
FGH,  the  inner  yield  curve  is  inoperative  and  plastic  flow  recommences 
at  H according  to  the  existing  flat  JCD. 

However,  if  the  unloading  trajectory  encounters  the  conical  portion 
of  the  yield  condition  as  at  P,  then  the  shrink  and  backtrack  technique 
must  be  used  to  create  the  new  yield  flat  RQU.  With  this  method,  although 
the  limiting  elastic  condition  during  elastic  behavior  may  be  non-convex, 
the  yield  condition  during  plastic  flow  will  always  be  convex.  For  example, 
regarding  Fig.  10  as  two-dimensional,  the  stress  trajectory  ABCEFG  produces 
the  non-convex  elastic  limit  consisting  of  straight  line  segments  TJ,  JE, 

EK,  KQ  and  the  circular  arc  QIT,  wheras  the  trajectory  ABCERS  produces  the 
convex  yield  condition  consisting  of  straight  lines  JE,  ER,  RU  and  the  circular 
arc  ULJ.  An  application  of  this  procedure  to  a plasticity  problem  with 
strain  hardening  will  be  presented  elsewhere. 
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[ Table  3:  Tip  displacements  and  stresses. 
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Table  4:  Plastic  strains. 


Table  6:  Unadjusted  plastic  strains 

and  yield  values. 
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Load-deformation  histories  for  f 
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Load  and  tip-deformation  trajectories 


Fig.  10 


